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Abstract 

It is difficult to predict how current climate change will affect wildlife species adapted 
to a tropical rainforest environment. Understanding how population dynamics 
fluctuated in such species throughout periods of past climatic change can provide 
insight into this issue. The drill (Mandrillus leucophaeus) is a large-bodied rainforest 
adapted mammal found in West Central Africa. In the middle of this endangered 
monkey's geographic range is Lake Barombi Mbo, which has a well-documented 
palynological record of environmental change that dates to the Late Pleistocene. 
We used a Bayesian coalescent-based framework to analyze 2,076 base pairs of 
mitochondrial DNA across wild drill populations to infer past changes in female 
effective population size since the Late Pleistocene. Our results suggest that the drill 
underwent a nearly 15-fold demographic collapse in female effective population 
size that was most prominent during the Mid Holocene (approximately 3-5 Ka). 
This time period coincides with a period of increased dryness and seasonality across 
Africa and a dramatic reduction in forest coverage at Lake Barombi Mbo. We believe 
that these changes in climate and forest coverage were the driving forces behind the 
drill population decline. Furthermore, the warm temperatures and increased aridity 
of the Mid Holocene are potentially analogous to current and future conditions 
faced by many tropical rainforest communities. In order to prevent future declines 
in population size in rainforest-adapted species such as the drill, large tracts of 
forest should be protected to both preserve habitat and prevent forest loss through 
aridification. 
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Introduction 

How a species has responded to climatic events not only in- 
forms us about the past but also offers insight into how the 
species might cope with future climate change (Davis et al. 
2005; Prost et al. 2010). The period during and since the Last 
Glacial Maximum (LGM) offers an ideal model for investi- 
gating this issue. The LGM, which occurred near the end of 
the Pleistocene at 26.5 Ka (Clark et al. 2009), was a dramatic 
period of glacial advance and global cooling and was fol- 
lowed by a period of warm and humid climate that lasted 
through the Early Holocene. Afterwards, fluctuations in 
aridity in the Mid and Late Holocene (along with human ac- 
tivity) shaped many of the Earth's modern ecosystems (Maley 
1996). Recently developed model-based analyses of molecu- 
lar data allow for the accurate inference of the timing and 
magnitude of past changes in the effective population size 
(N e ) of species through this time period (e.g., Beaumont 
1999; Drummond et al. 2005). The vast majority of studies 
that have performed these estimates have focused on taxa of 
either subtropical (Li et al. 2009; Liao et al. 2010), temperate 
(Depraz et al. 2008; Gratton et al. 2008; Lessa et al. 2010), 
alpine (Galbreath et al. 2009), arctic (Shapiro et al. 2004; 
Campos et al. 2010; Prost et al. 2010), or northern marine 
habitats (Larmuseau et al. 2009; Canino et al. 2010; Faurby 
et al. 2010; Marko et al. 2010). Very few studies have used 
model-based methods to investigate past changes in N e as- 
sociated with climate change in tropical species, particularly 
in relation to forest-adapted animals. This lack of knowledge 
is worrisome given the biodiversity present in tropical for- 
est ecosystems and the dramatic climatic and environmental 
changes that they are forecasted to undergo in the near future 
(the next 40-70 years; Beaumont et al. 2011; Heubes et al. 
2011). If we are to predict how climate change will affect 
species in forested tropical regions, we should have a better 
understanding of how they reacted to environmental changes 
during and since the Late Pleistocene. 

In an effort to address this knowledge gap, we investigated 
the past demography of a large-bodied terrestrial mammal 
adapted to a rainforest environment-the drill {Mandrillus 
leucophaeus) . The drill is an African Old World monkey (fam- 
ily Cercopithecidae) that is endangered, suffers from habitat 
loss and poaching, and ranks among the highest priorities 
for African primate conservation (Oates 1996). Like its sole 
living congenerer (the mandrill; M. sphinx), it is thought to 
have a relatively large home range and form multimale so- 
cial groups of up to 400+ individuals with females remaining 
in their natal group (Wild et al. 2005; Astaras 2009). Both 
species are also highly sexually dimorphic with adult males 
weighing over three times as much as females and possessing 
large canines and strikingly colorful perineal hair and skin 
(red, blue, violet) (Setchell et al. 2001). However, in contrast 
to the bright red and blue facial skin of male mandrills, male 
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Figure 1. Mandrillus leucophaeus male. 

drills have a jet black face with prominent cheek flanges and 
paranasal ridges framed by a rim of white skin and a red 
strip belowthe lower lip (Fig. 1). Also unlike mandrills, drills 
are primarily limited to a forested environment and are not 
known to exploit savanna-gallery forest mosaic ecosystems 
(Astaras 2009). 

The drill is endemic to the Cross-Sanaga-Bioko coastal 
forests of West Central Africa, which — as their name 
suggests — extend between the Cross River in Nigeria and the 
Sanaga River in Cameroon, and on Bioko Island (Equatorial 
Guinea) (Fig. 2). In the middle of the drill geographic range is 
Lake Barombi Mbo, which has one of the most complete pa- 
lynological records in the world since the Late Pleistocene for 
a forested region (Maley 1996; Maley and Brenac 1998; Bon- 
nefille 2007). Fossil pollen cores from this lake show a severe 
reduction in rainforest coverage between 24 and 1 1 Ka, which 
coincides with the LGM. This is followed by a dramatic forest 
expansion between 1 1 and 3 Ka when global temperatures, 
humidity, and percent forest cover were higher than current 
levels. At 3 Ka, the forest underwent another dramatic reduc- 
tion due to drier climate and increased seasonality. Between 
2 and 1 Ka the climate becomes more humid, which allows 
the forest to re-expand (Fig. 3; Bonnefille 2007). 
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Figure 2. Geographic distribution of the drill 
(Mandrillus leucophaeus) in the 
Cross-Sanaga-Bioko Coastal forests in West 
Central Africa (modified from Astaras 2009). 
Sampled localities are 1 = northern Cross River 
(several forests in the Boki region, including Afi 
Mountain and Okwangwo), 2 = southern Cross 
River (Western Oban Hills), 3 = Korup National 
Park, 4 = Ebo Forest, 5 = Bioko Island, Asterisk = 
Lake Barombi Mbo. 
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Figure 3. Fossil pollen profile for Lake Barombi Mbo displaying changes 
in forest coverage since the Late Pleistocene in the geographic range of 
Mandrillus leucophaeus (modified from Bonnefille 2007). Present day is 
on the left on the x -axis. Forest pollens decrease and savanna pollens 
increase during the Last Glacial Maximum (LGM) and during an arid 
period in the Mid-Late Holocene. Pollens from aquatic vegetation also 
rise during cold and/or dry periods as lake levels drop. 



Because it is restricted to a tropical rainforest habitat and 
exists in a region that has a documented record of environ- 
mental change during and since the LGM, the drill offers 
a rare opportunity to investigate how climate change might 
have affected a large-bodied tropical forest dwelling mam- 
mal. We hypothesize that the species underwent population 
size changes that mirror the changes in forest coverage in this 
region (based on the Barombi Mbo pollen profile); increas- 
ing and decreasing in population size as forests expanded 
and receded, respectively. In order to test this hypothesis, we 



sequenced ~2,000 base pairs of mitochondrial DNA from 
drills across the species' range and used a Bayesian Skyline 
Plot (BSP) to infer the timing and magnitude of past changes 
in drill female effective population size (N e {). 



Materials and Methods 
Sampling strategy 

We sampled 54 drills from Nigeria (northern and southern 
Cross River), Cameroon (Ebo Forest, Korup National Park 
[KNP]), and Equatorial Guinea (Bioko Island) between the 
years of 2005-2010 (see Fig. 2 and Table SI). The Nigerian 
localities included forests in the Boki region (e.g., Afi Moun- 
tain, Okwangwo) in the North and Western Oban Hills in the 
South. All together, the sampled localities encompass the far 
eastern, western, northern, and southern limits of the drill 
range, thus capturing the full extent of geographic variation 
in the species. Subpopulation sampling was also maximized 
to avoid incorrect inferences of population size change caused 
by recent immigrants from unsampled populations (Chikhi 
et al. 2010; Peter et al. 2010). 

All Nigerian and one of the Bioko samples were from wild 
caught captive individuals of known origin that are kept at the 
Pandrillus colony, while other samples were field collected. 
The vast majority of biomaterials were of fecal origin, al- 
though tissue biopsies from fingertips or liver were collected 
from some fresh carcasses during surveys of bushmeat off- 
take in the Malabo (Bioko) market in 2007. We specifically 
designed the sampling methods so they would not encour- 
age illegal hunting (i.e., no money was exchanged for market 
samples and groups were not habituated in the wild). 
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Genetic loci 

We analyzed two loci ("Brown" region and Cytochrome b) 
that together encompass 2,076 base pairs of the mitochon- 
drial genome. Mitochondrial DNA was chosen because it has 
a relatively fast mutation rate and a high copy number per 
cell, which makes it an ideal marker when using noninvasively 
collected samples that contain degraded DNA. The Brown re- 
gion is 935 base pairs long and includes the end of the ND4 
gene, beginning of the ND5 gene, and the intervening tRNAs 
(Brown et al. 1982). Cytochrome b (CYTB) is a 1,141 base 
pair long protein-coding region. Both of these loci have been 
used extensively in studies of Old World Monkey evolution- 
ary genetics and have been shown to provide an accurate 
account of female population history (e.g., Brown et al. 1982; 
Telfer et al. 2003; Newman et al. 2004; Wildman et al. 2004; 
Ting 2008; Zinner et al. 2009). 

DNA Extraction, PCR amplification, and DNA 
sequencing 

Total genomic DNA was extracted from fecal samples using 
the Qiagen DNA stool kit (Qiagen catolog number 51504) 
with modified protocols (overnight incubation in lysis buffer, 
half InhibitEX tablet, additional proteinase K, reduced elu- 
tion volume) in the labs of N. Ting and C. Roos. DNA was ex- 
tracted from tissue biopsies through a phenol chloroform ex- 
traction in the lab of N. Phillips soon after collection in 2007, 
stored at — 20°C, then processed through a whole-genome 
amplification in 2009 (Illustra GenomiPhi HY DNA Am- 
plification Kit; GE Healthcare Life Sciences). Amplification 
primers, sequencing primers, and cycling parameters for the 
Brown region and CYTB were from Newman et al. (2004) and 
Zinner et al. (2009), respectively. To ensure that these primers 
would not amplify nuclear pseudogenes in the drill, we com- 
pared Brown and CYTB sequences from one Bioko tissue 
sample that were sequenced from two different amplification 
protocols: (1) standard PCR and (2) a whole mitochondrial 
genome amplification via long-range PCR. This mitochon- 
drial genome was amplified in two 10,000 base pair overlap- 
ping halves, and the overlapping reads were sequenced and 
compared to ensure a circular molecule of mitochondrial ori- 
gin (following the methods of Thalmann et al. 2004; Raaum 
et al. 2005; primers available upon request). The Brown and 
CYTB sequences from the two different amplification proto- 
cols were identical, demonstrating that the standard short- 
range PCR primers are unlikely to amplify nuclear DNA. 
We also visually confirmed the coding frames and translated 
the protein-coding regions in our dataset in MacClade v4.08 
(Maddison and Maddison 2005) to check for unexpected stop 
codons, and we found no evidence that our primers ampli- 
fied nuclear pseudogenes. PCR products were cleaned using 
exonuclease I and shrimp alkaline phosphatase (Hanke and 
Wink 1994), cycle sequencing was performed using the Big 



Dye kit (Big Dye v3.1) following the manufacturer's proto- 
col for diluted reactions, and products were run on an ABI 
PRISM 3130xL or 3730 DNA Sequencer. Complementary 
strands were sequenced from multiple PCR products to en- 
sure the fidelity of the data, and the sequences were edited and 
assembled using Geneious v5.1.3 (Drummond et al. 2010). 
These data have been deposited in GenBank under accession 
numbers found in Table SI. Alignments of the Brown re- 
gion, CYTB, and a combined dataset were constructed using 
default values in ClustalW v2.0 (Larkin et al. 2007). 

Population structure and relationship 
analyses 

The presence of population structure can be mistaken for a 
genetic signature of population size decline. This is because 
many methods used to infer changes in past population size 
use models (e.g., Wright-Fisher, Moran) that do not account 
for structured populations. However, this is most problem- 
atic when there are cases of intermediate gene flow and recent 
immigration between subpopulations (Peter et al. 2010). In 
species with subpopulations that share very little gene flow, 
alleles coalesce in each subpopulation and a nonstructured 
population model may apply again (Chikhi et al. 2010). We 
expect this latter circumstance to apply to this study because 
we are using a matrilineal marker in a species that is likely 
female philopatric. With little movement of females between 
subpopulations, there should be a high amount of popula- 
tion structure among drill mitochondrial lineages. In order 
to assess this possibility, we used Arlequin v3.5 (Excoffier and 
Lischer 2010) to calculate Pst values among the Nigerian, 
KNP, Ebo Forest, and Bioko Island subpopulations using the 
combined dataset, conventional P-statistics, and 1000 per- 
mutations. In order to display the distribution of variation 
within the species, we constructed a haplotype network using 
statistical parsimony (Templeton et al. 1992) in the R package 
PEGAS (Paradis 2010). This method of network construc- 
tion has been shown to have lower error rates than many 
other network construction methods (Woolley et al. 2008). 
We also inferred a Bayesian gene tree and coalescent dates 
using the BEAST vl.6.1 package (Drummond and Rambaut 
2007). For this analysis, we conservatively discarded the first 
25% of trees from two BSP runs (see below) and combined 
them using LogCombiner. We then processed the tree file in 
TreeAnnotator and visualized the tree in FigTree vl.3.1. 

Evolutionary rate analysis 

We obtained separate Brown region and CYTB rates of evo- 
lution for use in the BSP analysis (see below) using both 
newly sequenced data and available nucleotide data in Gen- 
Bank from genera in the tribe Papionini. Brown and CYTB 
datasets were generated by aligning orthologous sequences in 
Papio hamadryas (Y18001), Theropithecus gelada (FI785426), 
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Lophocebus olbigena (JQ068156-Brown, JQ068153-CYTB), 
Cercocebus torquatus (JQ068155-Brown, JQ068152-CYTB), 
M. sphinx (JQ068154-Brown, JQ068151-CYTB), and M. 
leucophaeus (JQ068213-Brown, JQ068159-CYTB) using de- 
fault parameters in the program ClustalW v2.0 (Larkin et al. 
2007). The Lophocebus, Cercocebus, and Mandrillus sequences 
were previously unpublished data collected using the meth- 
ods of Raaum et al. (2005). Modeltest 3.8 (Posada and 
Crandall 1998) and PAUP* v4bl0 (Swofford 2003) were used 
to determine that the Hasegawa, Kishino, and Yano (HKY; 
Hasegawa et al. 1985) model was the best fit for both Brown 
and CYTB data under the Akaike Information Criterion ( AIC; 
Akaike 1974), with the former under a Gamma distribution 
(HKY + G) and the latter with a proportion of invariant sites 
(HKY + I) (Yang 1994). We used the BEAST vl.6.1 (Drum- 
mond and Rambaut 2007) package to estimate rates under an 
exponentially distributed relaxed uncorrelated clock model 
and a Yule speciation process tree model. Fossil data were 
used to set the divergence of Papio and Theropithecus to 5.0 ± 
0.5 million years ago under a normal distribution (95% CI: 
6-4 million years ago; Leakey 1993; Delson 2000). Four cat- 
egories of the Gamma distribution were used for the Brown 
region, and all other priors for model parameters and statis- 
tics were left at default values for both datasets. For each 
alignment, two 10,000,000 generation runs were conducted 
with parameters logged every 1000 generations, and log files 
were combined in LogCombiner after conservatively discard- 
ing the first 25% of samples as burn-in. Tracer vl.5 (Rambaut 
and Drummond 2009) was used to assess chain convergence 
and view the combined logfiles for both alignments, which 
yielded indistinguishable rates of evolution when rounded 
to the nearest 10 7 (3 x 10 -7 sub/site/generation; 9 x 10 -8 
lower 95% Highest Posterior Density; 5 x 10 -7 upper 95% 
Highest Posterior Density). Although female generation time 
(average age of reproduction; Kimura and Crow 1963) in wild 
drills is largely unknown, we estimate it to be 10-12 years. 
This number is derived using captive life history data from the 
Pandrillus colony in Nigeria (Wood 2007) on mean female 
age at first offspring birth (4.5 years), mean interbirth interval 
(1.4 years), and estimated lifespan (no more than 20 years, 
the age at which captive individuals deteriorate considerably; 
ELG personal observation). The average age at first birth and 
interbirth interval are similar to those from the semifree rang- 
ing mandrill (M. sphinx) colony at Lope (Setchell et al. 2001, 
2002, 2005), and a 12-year generation time is consistent with 
that which has been used in the closely related and similarly 
sized yellow baboon (Rogers and Kidd 1993). 

Effective population size change analyses 

We used DnaSP 5.0 (Librado and Rozas 2009) to calculate 
basic summary statistics on the combined Brown region and 
CYTB dataset. In order to assess past changes in female ef- 



fective population size (N e f ), we calculated Fu's F s (Fu 1997) 
and R2 (Ramos-Onsins and Rozas 2002). For both tests, we 
generated null distributions from 10,000 coalescent simula- 
tions of a constant-sized population to assess the significance 
of our results. We also used the BEAST 1.6.1 package to ana- 
lyze the drill Brown region and CYTB alignments through a 
BSP, which infers the timing and magnitude of past changes 
in population size (see Supporting Information for align- 
ments and the XML file). Because site and clock models were 
found to be the same in the within species alignments for 
both regions, we analyzed the two datasets under the fol- 
lowing linked substitution, clock, and tree models. The HKY 
model was found to fit both datasets best according to the 
AIC in Modeltest 3.8, and the clock model was set to a strict 
clock (because we believe rate heterogeneity in mitochon- 
drial protein-coding regions within a species such as the drill 
to be unlikely) with a 3.0 x 10 -7 sub/site/generation rate of 
evolution (12-year generation time). We also ran the analysis 
using a rate of 2.5 x 10 -7 sub/site/generation, which is the 
evolutionary rate if generation time is 10 years and average 
female drill lifespan is approximately 16 years. We used a 
Bayesian Skyline coalescent tree prior with 10 groups under a 
piecewise-constant model. Priors for model parameters and 
statistics were left at default values. The analysis was run for 
20 million generations with parameters logged every 1000 
generations, and Tracer 1.5 was used to inspect chain conver- 
gence and conduct the skyline reconstruction. We also ran the 
analysis twice to make sure results were consistent between 
runs and chains had converged. To assess whether the BSP 
generated the demographic scenario that best fits the data, we 
repeated this analysis under a constant population size tree 
model. A Bayes Factor Test (Kass and Raftery 1995) was then 
conducted in Tracer 1.5 on the likelihood traces to compare 
the fits of the BSP and constant population size models to the 
data (Suchard et al. 2001). 

We also conducted BSP analyses on numerous resampled 
datasets to determine the effects of unequal samples sizes 
among the four different populations (KNP n = 33, Nigeria 
n = 6, Ebo Forest n = 5, Bioko Island n = 10), assess potential 
population size changes in each population, and to investigate 
the effect of removing certain populations. The resampled 
datasets included KNP, Nigerian, and Ebo Forest populations 
individually (Bioko Island could not be analyzed alone using 
a BSP because it contained only one haplotype), Nigerian and 
Ebo Forest populations individually but randomly resampled 
to each containing 35 sequences, all populations but with 
the KNP population randomly resampled to 10 sequences, 
a KNP-Nigeria-Ebo Forest dataset (Bioko excluded), and a 
KNP-Nigeria dataset (Bioko and Ebo Forest excluded). All 
priors in the resampling analyses were kept consistent with 
the initially described BSP analysis. In addition, we repeated 
the original BSP analysis with an empty alignment to test the 
influence of the priors on the demographic results. 
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Table 1. Pairwise F ST values among sampled drill (Mandrillus leu- 
cophaeus) subpopulations computed using conventional F -statistics 
from haplotype frequencies (P < O.OB). 





Bioko 


Ebo 


Nigeria 


Korup 


Bioko 










Ebo 


0.845 








Nigeria 


0.884 


0.636 






Korup 


0.684 


0.556 


0.550 





Results 

Fsr values were all above 0.5 (P < 0.05; Table 1), indicat- 
ing that drill mitochondrial lineages are highly structured 
between Nigeria, KNP, Ebo Forest, and Bioko Island. The 
haplotype network (Fig. 4) shows one haplotype exclusively 
found in KNP, one haplotype shared between KNP and Nige- 
ria, a haplotype found exclusively in Nigeria that is closely 
related to a haplotype found exclusively in Ebo Forest, and a 




Figure 4. Statistical parsimony network of Mandrillus leucophaeus mi- 
tochondrial haplotypes. Five haplotypes are present in the sampled pop- 
ulations. Size of circle is proportional to frequency of each haplotype. 
Very little variation is shared between the different sampled populations. 



haplotype shared between Ebo Forest and Bioko Island. The 
gene tree (Fig. 5) shows close relationships between Ebo For- 
est and Bioko Island, and those two localities and Nigerian 
localities. All modern drill mitochondrial lineages coalesce to 
a date of 205,000 years ago (95% CIs: 162,500-331,000). 

Values of Fu's F s and R 2 were significantly different than 
expected under a model of constant population size (Table 2), 
with the observed values falling outside the 95% confidence 
intervals of the simulated values. The BSP indicates a nearly 
1 5-fold decrease in N e f that starts at the end of the Pleistocene 
and continues throughout the Holocene (Fig. 6). Specifically, 
median N et changes from 9,950 (95% CI = 4,112-28,777) 
to 681 (95% CI = 20-12,658). The decrease starts gradu- 
ally at the beginning of the LGM (~20-25 Ka), continues to 
the end of the Pleistocene and beginning of the Holocene, 
and becomes particularly sharp and pronounced in the Mid 
Holocene around 5 Ka. Between 2 and 1 Ka, the 95% CI 
becomes much larger and the mean N e f levels off slightly, 
although it still shows a declining trend. The Bayes Factor 
Test provides strong support for the Bayesian Skyline model 
when compared to a model of constant population size (log 10 
Bayes Factor = 3.131). The BSP analysis that used a shorter 
generation time showed a very similar result (Fig. SI). The 
resampled datasets all showed declines in Nef to varying de- 
grees except for the analysis of the KNP only dataset, whose 
mitochondrial lineages only allowed demographic inference 
of 10 generations into the past (Fig. S2). The analysis per- 
formed with the described priors and no data also resulted in 
a skyline with no change in N e {. 

Discussion 

Robustness of demographic analysis 

Our BSP analysis suggests a nearly 15 -fold decrease in drill 
female effective population size (N e f ) that began in the Late 
Pleistocene. We believe this demographic result to be statis- 
tically robust. The results of the F s and R 2 tests suggest that 
drill female effective population size (N e f ) may not have been 
constant in the past, and the results of the Bayes Factor Test 
demonstrates the BSP model was a better fit to the data than 
a model of constant population size. We also performed the 
BEAST analysis using the described priors with no data and 
recovered a flat skyline plot, thus indicating that the priors 
did not contribute to the inferred demographic history. 

While population structure can often be mistaken for a de- 
cline in population size, we do not believe this to be the case 
here. Signals of population structure and population decline 
are most commonly confused when there are intermediate 
levels of gene flow and divergent alleles are introduced from 
an unsampled subpopulation, thus mimicking the patterns 
of genetic diversity one expects from a population decline. To 
specifically avoid this, we employed a sampling strategy that 
incorporated the probable major subpopulations across the 
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Figure 5. Bayesian gene tree of Mandrillus leucophaeus mitochondrial lineages. Tree was inferred using a Bayesian Skyline Plot (BSP) tree prior in the 
BEAST v1 .6.1 package (Drummond and Rambaut 2007). Shown at major nodes are posterior probabilities, coalescent times (years, rounded to the 
thousand), and 95% confidence intervals. With the exception of two individuals (EB02, Glory), all individuals group with other members from their 
respective localities. 



Table 2. Summary statistics and results of tests for population size change as calculated in DnaSP S.O for the drill (Mandrillus leucophaeus). 
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number of polymorphic (segregating) sites; Cs 
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of synonymous polymorphic sites; Cn = number of nonsynonymous polymorphic sites; h = number of haplotypes; Hd — haplotype diversity; n = 
nucleotide diversity; jt 5 = silent nucleotide diversity; k = average number of nucleotide differences; F s = Fu's F s ; R 2 = Ramos-Onsins and Rozas R 2 
statistic; * = F s and R 2 results differed significantly from expectations based on the null model of constant population size. 



drill range (Chikhi et al. 2010). Furthermore, this issue is not 
as problematic when the population is extremely structured 
with little gene flow among demes (Chikhi et al. 2010). The 
very high F ST values among the drill subpopulations indicates 
this is the case in this dataset, as would be expected when using 
a matrilineal maker in a female philopatric species. In such 
a circumstance, it is likely that alleles coalesce in each sub- 
population and the inferred demographic collapse occurred 
across most demes. Our resampling analyses show that this 
was likely the case. Almost all resampled datasets resulted in 
a BSP that showed declines in N e { since the Late Pleistocene, 
with more severe declines occurring when sample sizes were 
higher. Removing certain populations or equalizing sample 
sizes among populations had little effect, and individual anal- 
yses of the Nigerian and Ebo Forest populations showed pop- 



ulation declines regardless of sample size. The lone exception 
is the KNP population, which failed to show a population 
size change. However, this result is not surprising because 
the mitochondrial lineages in this population had a recent 
coalescent time and only the past 120 years of demographic 
history were recovered in the BSP. 

Demographic history of the drill 

Because the study animal exists in a region with an excep- 
tional palynological record, past changes in N e can be directly 
correlated with known environmental change. Our hypoth- 
esis was that the drill would have experienced changes in 
population size that mirror changes in forest coverage. Our 
BSP analysis indeed suggests such a trend; a gradual decrease 
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Figure 6. Bayesian Skyline Plot (BSP) displaying changes in female effec- 
tive population size (A/ ef ) through time in Mandrillus leucophaeus based 
on 2,076 base pairs of mitochondrial DNA. Present day is on the left on 
the x -axis. A decline in M. leucophaeus A/ ef begins near the onset of the 
Last Glacial Maximum (LGM) and becomes most severe during the Mid 
Holocene (5Ka). 

in drill N e f that began prior to the LGM when global temper- 
atures were cooling, and the decrease persisted throughout 
the LGM when forests in West Central Africa saw a dramatic 
reduction. However, during the period of forest expansion 
following the LGM (10-3 Ka), we do not see evidence of a 
population size increase. JV e f continues to decline after the 
LGM and through the Early Holocene, and it abruptly de- 
clines in the Middle Holocene. There is, however, a slight 
leveling off of this decline in mean N e { between 2 and 1 Ka 
when the 95% CIs also become much larger. This last stage of 
the BSP with large margins for error may represent the limits 
of the mitochondrial data, where recent demographic events 
are difficult to detect. 

Causes of inferred population collapse 

Although changes to the environment can be a major in- 
fluence on population size dynamics, both competition and 
disease are known to affect animal abundance as well (Chap- 
man et al. 1999). However, we believe interspecific compe- 
tition was likely not a primary factor here because both the 
drill's large body size and social group size would provide 
significant advantages in competition for resources. While 
intraspecific competition (resulting from resource depletion) 
or disease remain possible causes, we believe environmental 
change induced by changes in climate to be the most likely 
primary explanation given that the timing of the population 
collapse coincides with known changes in forest coverage in 
this area. It should also be noted that these explanations are 
not mutually exclusive and can have a synergistic effect. For 



example, environmental disturbance can lead to stress, in- 
creased competition, and disease, all of which can act together 
to affect population dynamics in a species (Chapman et al. 
2006). 

Given that climate change is the likely primary driver of 
the inferred demographic history, it is surprising that there 
was no evidence of a drill population size increase during 
the period of warm, humid climate and forest expansion in 
the Early and Middle Holocene. This is counter to refugia 
theory expectations (Haffer 1982; Mayr and O'Hara 1986), 
which postulate that animal population range expansions 
and contractions should mirror climate-induced changes in 
habitat. This should particularly be the case in the drill, which 
has a highly diverse and adaptable diet and the ability to use 
resources in periods of drought that other sympatric primates 
cannot use, such as hard seeds (Astaras 2009; Astaras and 
Waltert 2010). A likely explanation is that our data do not 
provide sufficient resolution to detect a closer correlation 
between drill N e i and forest coverage. Use of a single locus 
may preclude inference of changes in N e that predate the most 
recent detectable demographic event, especially if that event 
was particularly severe (Heled and Drummond 2008). This 
is because a single locus will provide fewer coalescent points 
from which N e can be inferred compared to multiple loci. It 
is thus likely that the population decline seen in the BSP is 
primarily related to a reduction in forest coverage that began 
approximately 3 Ka (Fig. 3; Maley 2002). This general time 
period on the BSP is when the population collapse becomes 
most pronounced. Evidence of the collapse beginning earlier 
could thus be a residual effect of this event, a signature of an 
earlier demographic collapse related to the LGM, or a general 
lack of resolution. 

Whether or not the initial population decline was caused 
by the onset of the LGM, and whether or not there was a 
population expansion in the Early Holocene, it is clear that 
environmental conditions in the Mid-Late Holocene have had 
a dramatic effect on drill population dynamics and modern 
drill genetic diversity. Although some studies identified re- 
cent (i.e., the past few hundred years) habitat changes by 
humans as the main source of past population size change 
and current genetic variation in tropical mammals (Goossens 
et al. 2006; Olivieri et al. 2008; Craul et al. 2009; Thalmann 
et al. 2011), other studies are consistent with our conclusion 
that Mid-Late Holocene events were particularly influential 
(Heller et al. 2008; Okello et al. 2008; Lawler 2011). This 
time period saw an increase in aridity and the opening up of 
forests across Africa. Although this fragmentation facilitated 
the westward movement of Bantu-speaking populations into 
Central and West Africa, it is unlikely that human activities are 
responsible for the forest-cover changes reported for this area 
at the time (Maley 2002). Our results thus demonstrate the 
effects that climate change had on drill population dynamics 
around 3 Ka. 
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Implications for future climate change 

Modern increases in global temperature are already having an 
effect on wildlife communities (Hughes 2000; Walther et al. 
2002; Walther 2010), but predicting long-term responses to 
current climate change requires data from a multitude of 
sources. Investigating the effects of past climate change is 
thus one way of understanding how future climate change 
might affect the population dynamics and distributions of 
natural communities (Cordellier and Pfenninger 2009; Prost 
et al. 2010). The climate of the Mid-Late Holocene was char- 
acterized by an increase in aridity and seasonality during a 
warm interglacial period; and this is potentially analogous to 
current conditions where tropical climates are projected to 
experience increases in both temperature and aridity (Maley 
2002). This is due to global warming and human-mediated 
forest loss that disrupts the moisture cycle and may lead 
to increased drought and seasonality (Hannah et al. 2007; 
Hannah and Lovejoy 2007; Heubes et al. 201 1). Based on the 
findings of this study, forest-adapted species such as the drill 
may not be able to cope with such changes and their popu- 
lations could dramatically decrease — a catastrophic scenario 
for a species already on the brink of extinction. Limiting for- 
est loss by protecting large tracts of intact tropical forest will 
therefore pay dividends beyond just the protection of wildlife 
habitat from poaching and human development — it will aid 
in maintaining the evapotranspiration cycle and potentially 
prevent an increase in aridity and seasonality that could have 
negative effects on species such as the drill. 

Conclusions 

We conducted a study of changes in past effective population 
size since the LGM in a tropical rainforest-adapted mammal 
to see how it responded to periods of past climate change. 
Our study animal was the drill, which is a poorly known 
and endangered primate (Oates 1996), and this research is 
the first genetic work to be conducted on this species. We 
predicted that drill population dynamics would closely mir- 
ror changes in West Central African forest coverage since the 
Late Pleistocene. Our data suggest that the drill underwent 
a population collapse that was most severe during the Mid 
Holocene, which supports the conclusions of other studies 
of sub-Saharan mammals that found events during this pe- 
riod to have had a dramatic impact on population dynamics 
and modern genetic diversity (Heller et al. 2008; Okello et al. 
2008). In West Central Africa, this coincides with a reduc- 
tion of forest coverage associated with warm temperatures 
and increased aridity and seasonality, and is potentially anal- 
ogous to the foreseen global warming effects in the region. 
In order to prevent or at least curb future population re- 
ductions in rainforest-adapted species such as the drill, large 
areas of forest must be protected to both preserve natural 
habitat and prevent disruptions in the moisture cycle. Future 



research that incorporates more genetic loci and multiple 
sympatric taxa will provide better resolution in reconstruc- 
tions of demographic history and shed further light on how 
species might respond to the same ecological changes. Such 
information is essential for the better design of landscape 
scale — rather than species focused — conservation strategies 
in the region. 
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Supporting information 

Additional Supporting Information may be found online on 
Wiley Online Library. 

Figure SI. Bayesian Skyline Plot displaying changes in fe- 
male effective population size (N e f) through time in Man- 
drillus leucophaeus based on 2,076 base pairs of mitochon- 
drial DNA, a 10 year generation time, and a 2.5 x 10~ 7 
sub/site/generation mutation rate. Present day is on the left 
on the x-axis. A decline in M. leucophaeus N e f begins near the 
onset of the Last Glacial Maximum and becomes most severe 
during the Mid Holocene (5Ka). 

Figure S2. Bayesian skyline plots of resampled data. Shown 
are the median and 95% confidence intervals. Time is on x- 
axis with units in generations and a maximum value of 5000, 
except for d) which has a maximum value of 20. Female 
effective population size on y-axis with a maximum value of 
100,000 individuals, a) All populations with Korup National 
Park (KNP) population resampled to 10 individuals, b) KNP, 
Nigeria, and Ebo Forest populations (Bioko Island excluded), 
c) KNP and Nigeria populations (Ebo Forest and Bioko Island 
excluded), d) KNP population only (generation time on x- 
axis with a maximum value of 20). e) Ebo Forest population 
only, f) Ebo Forest population resampled to 35 individuals, 
g) Nigeria population only, h) Nigeria population resampled 
to 35 individuals. 

Table SI. Sampled individuals with associated locality data 
and Genbank accession numbers. Biomaterials from Bioko 
individuals with an "Mleu#" ID were fingertips (115, 120, 
135, 160, 37) or liver (194, 203, 205) from fresh bushmeat 
carcasses. All other biomaterials were of fecal origin. 

Please note: Wiley- Blackwell is not responsible for the content 
or functionality of any supporting materials supplied by the 
authors. Any queries (other than missing material) should be 
directed to the corresponding author for the article. 
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